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ABSTRACT 

The body of work presented here revolves around the investigation of the existence 
and nature of extra-solar planetary systems. The fitting of stellar radial velocity time 
series data is attempted by constructing a model to quantify the orbital properties 
of a star-planetary system. This is achieved with the Planetary Orbit Fitting Process 
(POFP). Though specific to the investigated problem, the POFP is founded on two 
separate, more general ideas. One is a Solver producing the gravitational dynamics of 
a Three-Body system by integrating its Newtonian equations of motion. The other is 
an independent optimisation scheme. Both have been devised using MATLAB. Ap- 
plying the optimisation to the Solver results in a realistic Three-Body dynamics that 
best describes the radial velocity data under the model-specific orbital-observational 
constraints. Combining these aspects also allows for the study of dynamical instability 
derived from interaction, which is reaffirmed as a necessary criterion for evaluating the 
fit. The validity of POFP solutions with respect to the observations and other mod- 
els is discussed in this context. The underlying generality and fundamental principles 
demonstrate a larger frame of operation where problems in Physics and Mathematics 
can be solved with a multitude of techniques. 
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INTRODUCTION 



One way to infer the existence of objects orbiting a star 
is by the use of the stellar radial velocity (RV) measure- 
ments obtained through the observation of the Doppler Shift 
in the star's spectrum. It can be shown (Windmiller 2006) 
that in a n-Body gravitational system in which the centre 
of mass (CM) does not accelerate, the secondary centre of 
mass of n — 1 bodies tends to move with a momentum in 
the opposite direction to that of the other body, thus re- 
sulting in a reflex motion that has a component along the 
observer's line of sight. It is therefore possible to construct 
a model in which the stellar motion in the radial direction 
is calculated based on the gravitational dynamical interac- 
tion of the bodies in such a system. Numerous stars have 
been surveyed to date by different observing groups (UCLES 
Keck/Lick, HARPS /CORALIE, The AAPS...) using high 
resolution velocity measurements of each star with plane- 
tary bodies either determined or conjectured. The method 
of analysis of these radial velocity data has been outlined in 
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several papers over the years (Laughlin 2001, Rivera 2005, 
Pepe 2006) and has slowly evolved with the essential fea- 
tures remaining unchanged. The first step in the analysis 
is the determination of Keplerian orbital parameters under 
the assumption of no mutual gravitational perturbations be- 
tween the planets, eliminating any significant periodicities in 
the residual. Once the number of conjectured planets - each 
with its own Keplerian parameters - has been determined, 
additional assumptions such as co-planarity, edge-on orbits, 
or resonance, are used to arrive at a full set of Newtonian 
initial conditions. This orbital configuration is then used as 
the initial guess for a quasi-Newton's method search (usu- 
ally the Levenberg-Marquardt algorithm). The restriction of 
no mutual planetary gravitational interaction is lifted, and 
the search proceeds to a self consistent local minimum in the 
parameter-space vicinity of the Keplerian orbital fit. Using 
this family of methods, only a local minimum close to the 
Keplerian parameters under the specific restrictions will be 
found. 

The alternative model proposed here for the prediction 
of stellar RV data relies on a direct integration of the New- 
tonian equations of motion for a Three-Body gravitational 
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interaction. The developed RV reduction algorithm consti- 
tuting this model encorporates the assumption that in gen- 
eral, there will be significant mutual gravitational planetary 
perturbation, since it is exactly the star's reflex action to the 
planetary presence that is being measured. Searching the 
full set of possible initial conditions while directly integrat- 
ing allows for the investigation of all orbital pathologies un- 
constrained by any factors (such as co-planarity, resonance, 
etc.). Each solution, which is the set of initial conditions for 
the single integration of Three-Body trajectories and veloc- 
ities is put in the form of a 12-tuple of orbital parameters. 
Once the full Newtonian trajectories are calculated, an algo- 
rithm based on the Singular Value Decomposition finds the 
optimal orientation based on the RV data set. In order to 
find a solution which predicts a RV curve that sufficiently 
matches the observations, an optimization process based on 
methods from the literature has been devised. The search for 
the absolute minimum of the \ i s implemented first with 
a Genetic Algorithm, then with a local search method, typ- 
ically a Compass or a quasi-Newton's method. The entire 
process of analyzing and predicting the RV data, is referred 
to here as the Planetary Orbit Fitting Process (POFP). 
All of the algorithms described here are written in MAT- 
LAB, utilizing its versatility as a mathematical computer 
language, and a working platform with a graphical interface 
and powerful application packages. These include diverse 
tools such as the suite designed to integrate and manipulate 
Ordinary Differential Equations (ode). Because the POFP 
in its entirety was not previously described in the literature, 
a comprehnsive description is included here, complimented 
by the appendices for more detail. The Three-Body Solver 
at the centre of the POFP is described in §2. The application 
of the Solver to produce a predicted stellar RV time series 
via the orbital-observational parameters is described in §3. 
A brief description of the optimization scheme follow in §4. 
Considerations involving the comparison between predicted 
and observed quantities are discussed in §4 as well. The Ke- 
plerian model which is the basis of the methods currently 
used throughout the literature, orbital stability considera- 
tions and the connection between the two are discussed in 
§5. Applications of the POFP using the published radial ve- 
locity data for the HD160691 and HD141929 systems are 
demonstrated in §6 and other cases are discussed in future 
papers. 



2 THE SOLVER 

The gravitational interaction between N bodies due to New- 
ton's Second Law and the Universal Law of Gravitation is 
given by the system of equations 
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where rrij is the mass of particle j applying force on particle 
i from a relative position nj. In the case of three bodies in 
three dimensions, N = 3, and ri — (x, y, z) is the displace- 
ment vector of particle i. Assuming the CM of the system 
moves with constant velocity, one of the vectors in system 
(1) can be expressed in terms of the other two: 
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where r is the acceleration vector in the CM frame. The 
system (1) then consists of six second order, or twelve first 
order differential equations of motion. The initial conditions 
for this system are the three masses and their initial position 
and velocity vectors. 

It is then integrated for the trajectories ri(ti), and ve- 
locities Vi(U) utilizing the MATLAB ODE suite. The al- 
gorithms chosen from the suite were odell3 (Variable or- 
der Adams-Bashforth-Moulton) and ode 45 (Runge-Kutta- 
Fehlberg, Dormand-Prince pair). The time span for integra- 
tion, the relative and the absolute tolerance were specified 
along with the initial conditions. As part of this scheme, 
interpolation can be made over the integrated interval to 
any desired step siz^B allowing for no further integration. 
The Solver has been tested on problems with known so- 
lutions to make sure of its proper use. With the integra- 
tion of every orbit, the angular momentum is calculated for 
each time step and the relative error in it throughout the 
entire integration time is computed to show conservation. 
Using the constants defined to great accuracy by the IAtQ, 
and data tabulated in Nautical Almanac Office (2006), the 
trajectories of the Earth-Moon Barycentre, the Earth-and- 
Moon and the Earth-and- Jupiter systems around the Sun 
were computed. The orbits were shown to close (in distance 
and time) to a sufficient accuracy compared with that of 
the initial data (Windmiller 2006). All orbits showed the 
conservation of angular momentum to the number of digits 
defined by the error-tolerance. The trajectories calculated 
by the Solver have no specific orientation towards an ob- 
server. To define the direction of the observer's line of sight, 
the Solver's coordinate system must be rotated (see Ap- 
pendices B,C). 



3 APPLICATION OF THE SOLVER TO A 
PLANETARY ORBITAL MODEL 

The application of the Three-Body solutions (trajectories 
and velocities as a function of time) produced by the Solver 
to the problem of fitting an observed stellar RV variation, is 
made through the assumption that two of the masses (plan- 
ets) are orbiting the third (star). The search of the full set 
of initial conditions for the Newtonian equations of motion 
must bring into account the fact that any random set of such 
conditions most likely leads to unbound trajectories. To re- 
solve this problem, the set of feasible initial conditions can 
be parametrized by the set of instantaneous Keplerian pa- 
rameters at the time of the initial data point. Masses, eccen- 
tricities and angular information as to relative positions of 
the orbits and their periastra are used for this parametriza- 
tion. The result is a set of orbital parameters from which 
the initial conditions for the orbiting masses in their two 
graviataionally bound trajectories are calculated. 

In the plane of each orbit, the initial position and ve- 
locity vectors are defined with respect to the point of pe- 
riastron, using the eccentricity, e, and a 'deZaj/' angle, 0, 
measured from the semi-major axis (see Appendices A,E). 
Given the relative positions of the two semi-major axes and 

1 Details on the workings of the ode suite can be found in 
Shampine et al. (2003) 

2 http://193.49.4.189/Units.234.0.html 
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the relative tilt between the orbital planes, the initial posi- 
tions and velocities of one orbit can be transformed into the 
coordinate system of the other so that all quantities are de- 
fined in the same system. This transformation is achieved by 
a multiplication of three rotation matrices in which the an- 
gles are defined by the relative coordinates of the semi major 
axes (see Appendix A) . These nine orbital parameters allow 
for any possible geometry of the two planets orbiting the star 
at the epoch of initial conditions. The observer's coordinate 
system is defined as one where the (x' , y') plane is the plane 
of the sky and the z' is the radial direction. The orientation 
of the Solver's coordinate system to the observer's line of 
sight can be defined using two angles. One of the orbits is 
chosen and its inclination, i, and longitude of periastron, ui 
are used to transform the coordinate system of the Solver 
into that of the observer. This transformation is also the 
product of two rotation matrices (Hilditch 2001). The orbital 
parameters that are therefore sufficient to define the initial 
conditions to the Solver for the planetary model are the ec- 
centricities ei^, the delay angles #1,2, the length of the first 
semi-major axis x\ and the coordinates of the second semi- 
major axis (2:2,3/2,22), the relative orbital tilt angle r, and 
the three masses mi, 7712,7713. The observational parameters 
required are the inclination i and the longitude of perias- 
tron u). Parameter-space for any n-tuple of initial conditions 
satisfying the Three-Body requirements of the Solver has 
fifteen dimensions. This 14-tuple of planetary-model param- 
eters thus reduces by one the number of dimensions in which 
the search for initial conditions must take place. 

The number of dimensions can be reduced further to 
include the 12 orbital parameters only. For any orbital con- 
figuration defined by the parametric 12-tuple, the observer's 
position (i.e. i and ui) from which the Solver-frame stellar 
velocity best describes the RV data can be calculated. The 
calculation is based on optimizing a least square solution 
for a unit 'look vector', defined by the two angles, which 
points from the CM of the stellar system to the observer. 
The matrix equation is solved using the singular value de- 
composition (utilizing MATLAB's svd function) applied on 
the star's velocity in the Solver frame (see appendix C). 



4 OPTIMISATION AND THE FITNESS 
FUNCTION 

To find the parameter-point (14-tuple or 12-tuple) for which 
the comparison between the model-predicted and the obser- 
vational data results in the absolute minimum y 2 , an op- 
timization process was carried out in two stepqj. The first 
step is the implementation of the Genetic Algorithm (GA) 
as a global search for the region in parameter space where 
the lowest x 2 ues - The second step of refining the search 
to the exact parameter-space location of the absolute min- 
imum in the y 2 was carried out in various ways. Of the 
derivative-based methods, the linear slope and the quadratic 
slope methods were found to be the least efficient. This is 
because of the rapid change in the fitness function's slope 
and numerical-derivative noise problems ( Windmiller 2006) . 
The Hybrid Newton-line search method was typically used 

3 All algorithms discussed here were developed by Donald Short 



to advance the search quickly where other methods were un- 
satisfactory. The second-step method most frequently used 
was the Compass search. All algorithms were constructed 
as independent MATLAB programs. The programming of 
the GA followed several discussions in the literatur^f|[f|. The 
compass was based on Kolda, Lewis & Torczon (2003). All 
methods used, most specifically the GA and Compass, were 
tested against known analytical functions before their use in 
the POFP. 

To produce the predicted radial velocity data, the 
POFP implements two translation processes. First, the n- 
tuple of parameters is translated into initial conditions for 
the Solver. The information about the dynamical proper- 
ties of the system (in the form of the solution produced by 
the Solver) then undergoes a translation into the observer's 
plane of the sky coordinate system. An n-tuple of parame- 
ters therefore constitutes a single hypothesis that a specific 
configuration of initial conditions results in a model that is 
compatible with the observations. The quantification of this 
level of compatibility is made via a fitness function. The 
specific fitness function used by the POFP is the y 2 (Chi- 
Squared). The derivation and treatment of the y 2 function 
follows that of Bevington & Robinson (2003) and Press et al. 
(1992). The assumption that the CM of the stellar system 
does not accelerate allows for a constant offset velocity to 
be incorporated into the y 2 function (such a fitness function 
is presented in Goidziewski & Maciejewski (2003)). The op- 
timal offset velocity is calculated in closed form within each 
evaluation of the y 2 (appendix D). Because of technical con- 
siderations and as the standard used in the literature, the 
form of the fitness function that is used in practice is the 
reduced Chi Squared, y 2 > (with v being the number of data 
points minus the number of parameters, or degrees of free- 
dom, plus 1). The expected value of y 2 is ~ 1- However, 
because this exact value typically is not achieved, a range 
of values must be considered. One approach in the interpre- 
tation of the POFP solutions' y 2 is based on using the F 
test as suggested in Press et al. (1992). Given xt x an d y 2 2 it 
is not possible to reject the assumption that both represent 
the same data set if both reside inside the same interval. 
This interval is determined as a confidence measured by a 
probability which is given by the appropriate integral-test 
probability table. For a desired probability of not exceeding 
the ratio F = xtt /y 2 2 (typically a probability measuring the 
confidence of one standard deviation), the ratio can be cal- 
culated, considering the numbers of degrees of freedom ^1,2- 
Press et al. (1992) provides the expression for calculating the 
ratios as a function of the probabilities. This interpretation 
allows for both a comparison between any two values of y 2 
(provided v is similar in both and sufficiently large) and a 
determination of the range of values in which solutions are 
acceptable relative to one another. 

When considering the acceptability of a given solution 
produced by the POFP, the xt is treated as the primary 
criterion. The secondary measure for how well the fit de- 
scribes the observations is the Root Mean Square (RMS) of 



4 Houk, Joines and Kay: http://www.ie.ncsu.edu/mirage/ 
G AToolBox / gaot /papers/ gaotv5.ps 

5 Parker: http://pharos.cpsc. ucalgary.ca/Dienst/Repository/2.0/ 
Body/ncstrl.ucalgary_cs/1995-564-16/pdf 
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the error, along with the distribution of the residual veloci- 
ties which is expected to be a Gaussian centred at zero. In 
the framework of the POFP, a useful RMS-related quantity 
is RMSerr = -\/af cs — <r? it . Here, a ICB is the standard devia- 
tion in the residual velocities (comparable to the RMS in the 
literature) and ay lt is an uncertainty added to the observed 
velocities due to a variety of factors largely unknown and, 
in part, related to stellar chromospheric activity (Wright 
2005). The latter is termed here Jitter. Apart from consid- 
ering its magnitude, the quantities RMS and RMSerr can 
be compared to the uncertainty in the observed data. The 
quantification of the fit by the xt an d the RMS may yet be 
insufficient for determining how well the model describes the 
observed phenomenon. One of the most important consider- 
ations that typically applies is that of stability. As the inves- 
tigation using the POFP of the 47UMa system (described in 
Paper II), as well as that of HD160691 show, solutions that 
have a low xt may also ascribe a large eccentricity to one 
of the planets, resulting in instability. Dynamical stability 
of the system, therefore, must be considered as a tertiary 
criterion for accepting the fit. 



5 THE KEPLERIAN MODEL AND STABILITY 

The application of the POFP is demonstrated in §6 on the 
star systems HD160691 and HD141929 which have previ- 
ously been investigated in the literature. In order to repro- 
duce published results that are based on a Keplerian model, 
a similar scheme has been developed as part of the POFP 
programming. Separate from the Planetary Model described 
earlier, this Keplerian Constants (KC) model allows for a 
confirmation or an evaluation of the published Keplerian re- 
sults and a comparison between those and the ones produced 
by the POFP. 



5.1 The Keplerian Constants Model 

The classical Keplerian model that has been programmed 
here is termed the KC model. It is based on the equations 
derived for binary mechanics as found in any standard text- 
book (e.g. Hilditch 2001) in relation to the general analytic 
Two-Body solution and is comparable to that commonly 
used in the literature. In this model, the motion of a mass 
in an elliptical orbit around a second mass (here a planet 
around a star) as a function of time is described by Ke- 
pler's equation. Using the period, P, and the periastron 
passage time, T, this equation relates the Mean Anomaly 
O = 2n{t — T)/P of the orbiting mass to its Eccentric 
Anomaly, E, and the eccentricity of orbit, e: 



V r (t) = K [cos (0(t) +lo) + ecosw] 



(5) 



9 = E - e sin E 



(3) 



The program iterates (using the Newton Raphson method) 
for the Eccentric Anomaly. The true anomaly (delay angle) , 
9, is calculated as a function of time: 



9{t) = 2arctan 



1 + e 



■ tan 



E(t) 



(4) 



Using the Half Amplitude of the observed velocity, K, and 
the Longitude of Periastron of orbit, w, the radial velocity 
can be predicted: 



The predicted radial velocity is calculated separately, using 
Equation (5), for the influence of each planet on the star. 
These contributions to V r are then superimposed to give 
the total predicted radial velocity. The planets are there- 
fore assumed to have no interaction, making the Keplerian 
quantities constant in time. A parameter-point constituting 
a solution in this model is the n-tuple of KC (K, P, e, T, w) ■ 
with j designating the planet (i.e. the 5-tuple j = 1 for a sin- 
gle Keplerian, the 10-tuple j = 1,2 for a double Keplerian, 
etc.). The appropriate parameter-point which minimizes the 
xt can be searched for using the optimization scheme de- 
scribed in §4. As with the POFP, the GA is typically used 
as a first step, followed by the Compass method. When opti- 
mizing, the KC model is highly sensitive to the initial guess 
(or initial population). 

POFP integration over a sufficiently long time span 
(typically much longer than that of the period of the outer 
planet) can expose a system as being dynamically unstable 
or prone to instability. This is done directly by showing an 
ejection of one of the bodies from the system, or at the very 
least a strong gravitational influence of the two planets on 
one another, making noticeable changes to their trajectories. 
Because the KC model assumes no gravitational interaction 
outside that between the star and each planet, a separate 
analysis for stability is required where other factors, not in- 
herent to the model, are brought into account. An example 
of a discussion involving this type of analysis can be found 
in Goidziewski & Maciejewski (2003). Using the Solver to 
compute the trajectories makes for a direct integration of 
the Newtonian equations of motion which is general in its 
underlying assumptions, and geometry. Because of that, a 
translation of the initial conditions of the KC model into 
those of the POFP, followed by integration based on the 
latter, may be useful in approaching the problem of orbital 
stability. 



5.2 From Keplerian to POFP 

The KC parameter point (K, P, e, T, uj)- used by the model 
to express the orbital properties of the j th planet as they 
appear to an observer through the plane of the sky (i.e. in 
the form of stellar RV data), describes a family of Newto- 
nian solutions, each being an orbit initiated by a different 
set of position and velocity vectors. These solutions differ 
from one another by two parameters. The first is the plan- 
etary mass traveling in each orbit which can alternatively 
be expressed by the inclination of that orbit (using P, K, e 
and the mass of the star). For any two orbits, the second 
quantity is the difference in their position angle of the nodal 
point, (usually designated by fi). Therefore, the Newtonian 
fit for any Keplerian solution can be described by the xt 
calculated as a function of the masses or inclinations of the 
planets and their differences in position angles. Because the 
POFP's primary use is in solving for three-body trajecto- 
ries, a rigorous analysis of the special case of two planets 
lends itself as a first step in understanding the relationship 
between the different models. Given the stellar mass, a 10- 
tuple of double Keplerian parameters (i.e. 10 KC) can be 
translated into the POFP 12-tuple of parameters, excluding 
the planetary masses. The xt for such a solution can then be 
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calculated for triplets of (11,12, Af2), i.e. for pairs of masses 
(mi, 7712) and their orbits' relative positions in the plane of 
the sky. As a preliminary analysis of this three dimensional 
parameter space, a crude optimization can be made where 
the AO, that minimizes the xt f° r each pair of masses is cho- 
sen. The resulting surface, xt = s ( m ii m 2) An . can then 
be examined. For the entire span in inclination, the domain 
of this surface shows the allowable range of masses the par- 
ticular Keplerian solution can have. The regions of lower xt 
may easily be located on this surface, and the stability of 
solutions represented by specific points in them can be in- 
vestigated. The principal idea of such an analysis through a 
translation is shown here in §6 for the double Keplerian as 
it applies to the test case of HD 160691. Some detail of the 
translation is given in appendix E. A more rigorous stability 
analysis based on these methods will be discussed in a paper 
to follow. 

5.3 Investigating Stability 

One qualitative argument regarding interaction in a stellar- 
planetary system is the following: if a gravitational influence 
of the outer planet on the star exists (one that is sufficiently 
strong to be detected) then such an influence of the outer 
planet on the inner planet must also occur, causing the plan- 
ets to interact, provided that the two are on the same side of 
the star. Concretely, in order to assess whether the trajec- 
tory calculated from a set of initial conditions results in an 
ejection of one of the bodies, or more generally, is prone to 
instability, the Solver can be used to integrate the orbit on 
a sufficiently large time span. Such an evaluation, based on 
the graphical representation of the trajectories, was carried 
out for all of the POFP solutions in Table 1 and Table 2 dis- 
cussed in §6. The typical time span used was on the order of 
hundreds of times the larger period of the two planets. The 
assessment, using this method, of a possible instability of 
an orbital configuration requires a large degree of intuition 
on the one hand, and is limited by computational resources 
on the other hand. Because of that, the Hill stability con- 
dition was applied as a complementary criterion, following 
Gladman (1993). In his paper, Gladman (1993) distinguishes 
the Hill stability, where the bodies are forbidden from un- 
dergoing close approach for all time, from the stronger La- 
grange stability where, in addition, the semi-major axes are 
bounded. The author defines closest approach, which is re- 
garded by him as instability, as being inside the sphere of 
influence of the outer planet where "...it is better to view 
the motion as one planet around the other being perturbed 
by the 'Sun' rather than the planets moving around the Sun 
and being perturbed by each other" (Gladman 1993). The 
cutoff for this closest approach is defined numerically by the 
author in terms of the planetary mass ratio. As is the case 
in Gladman (1993), the approach taken here for establish- 
ing the Hill instability relies on the parametrization of the 
Three- Body dynamics by Marchal & Bozis (1982) in terms 
of the three masses, the gravitational constant, the angular 
momentum and the total energy of the system. For the case 
where two of the bodies are much smaller in mass than the 
third one, the critical value of this parametric expression, in- 
dicating a bifurcation of the system's dynamics, is computed 
by Marchal & Bozis (1982) in terms of the masses alone. The 
difference between the Hill number of a particular solution 



(i.e. the dynamical parameter of that orbital configuration) 
and the Hill critical value is referred to here as the Hill dif- 
ference. When the Hill difference is positive, the planets are 
forbidden from having a disruptive close approach, which 
means the system is Hill stable (although ejection of the 
outer planet may still occur). Fulfilling the Hill condition is 
stated by Gladman (1993) to be sufficient but not necessary 
for overall stability in the system's dynamics. 



6 DEMONSTRATING THE POFP 
APPLICATION 

The case of the double-lined binary HD141929, observed and 
investigated by Carrier (2002), was used to show how the 
POFP treats a degenerate Three-Body structure. The anal- 
ysis of this stellar system included a reproduction of the 
published results using the KC model, at least two viable 
Newtonian solutions competing with the published and re- 
produced ones, an error analysis on all solutions, and a xt 
surface showing the expected stability properties of the sys- 
tem. Because the process was started with as few assump- 
tions as was possible about the system (such as the stellar 
mass of the primary derived from its spectral type), a rela- 
tively general initial population was given to the GA. The 
solutions produced were all in the form of a translating bi- 
nary with a small third mass at a distance, typically ejected 
from the hierarchy. These solutions can then be optimised 
further with the third mass set to an insignificant value. The 
error analysis of the KC reproduction, showing values close 
to the published ones is showin in Table 3. Some of the main 
results for this system are shown in Table 4. 

6.1 An application to the star HD160691 

To demonstrate the intricacies of the POFP analysis and 
its intepretation, the more complex case of HD 160691 was 
chosen. Although more RV data of this system have accumu- 
lated since the publication referred to in this paper, the case 
of HD160691 is brought as an example of an analysis which 
is both comperable to other publications at a given point of 
time, and a work in progress. Some of the conclusions emi- 
nating from this analysis, especially in light of the new data, 
allow for the exploration of new directions in research. 

The RV observations of HD160691 and a double Ke- 
plerian fit to these data are provided by McCarthy et al. 
(2004). The authors discuss stability analyses done by Bois 
et al. (2003) and Goidziewski & Maciejewski (2003). Tables 
1 and 2 here show the published result, its Keplerian re- 
production, and the POFP solutions for this star, in a form 
compatible with the literature. For all solutions generated 
here, it is assumed that the data were taken using a single 
instrument. Table 1 contains the stellar mass and the quan- 
tities for the first planet, along with the xi an d RMS of 
each filQ, and the stellar offset velocities. Table 2 contains 
the properties of the second planet and comments regard- 
ing both planets. In all of these solutions, the stellar mass 

6 The RMS for the solutions presented here is the standard devi- 
ation of the residual, (<r ros ) without subtracting the uncertainty 
due to jitter (<Tjit) from it, making it higher than RMSerr dis- 
cussed in §4. 
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was not allowed to vary by more than 15 percents from the 
published value of 1.08Mq. The maximum radial velocity 
for each planetary mass m p around the star was calculated 

as 

_ GnipP sin i 

Solution #1 is that published by McCarthy et. al. 
(2004). Using the solution's published parameters in the KC 
model described in §5.1, without including any jitter in the 
uncertainties, results in xl = 385 with RMS = 47.7. So- 
lution #2 is a reproduction of solution #1 using the KC 
model and applying the optimization described in §4. The 
published parameters (yielding solution #1) are used as a 
basis for an initial population and the GA followed by the 
Compass are implemented to find the minimum in the fit- 
ness function, with a result of xl = 1-27 with 3 m/s jitter, 
or xl = 2-7 with the jitter removed, and RMS = 4.2. The 
high xl m solution #1 could be attributed to an error in one 
or more of the published parameters, most specifically the 
published value of the periastron passage of the first planet. 
Substituting this value with the guess, T = 260, based on 
the optimization in solution #2, yields a new xl = 1-26 with 
3 m/s jitter or xt = 2.68 (xv = 1.63) without any jitter, and 
RMS = 4.15. The published solution is therefore adequately 
reproduced. 

Solution #3 is the first one listed that was generated by 
the POFR With a xl = 1-44 and a RMS = 4.30, it shows an 
agreement with the McCarthy et al. (2004) parameters for 
the first planet. Comparing the fit of solution #3 to that of 
the reproduction solution #2 (both include the 3m/s jitter 
in the data) shows the xl °f #3 to be larger than that of 
#2. However, as explained in §4 and is shown in Table 5, the 
two xl values are well within a la interval of one another. 
The RMS of #3, a measure of the fit which does not depend 
on the size of the error bars, is also larger than that of #2. 
The small negative Hill difference for this case could indicate 
instability. The solution was integrated over a time span of 
400 times the larger planetary period. Examining the plot 
shows the trajectories to have a possible instability, i.e. not 
a definite ejection of one of the objects, yet a relatively large 
degree of interaction that may be disruptive subsequently. 
The period of the second planet in #3 is larger than that 
of #2 and the eccentricity much smaller. Nevertheless, it is 
possible that the large mass of the first planet, M P1 = 63Mj 
is what makes this solution prone to instability. 

Solution #4 demonstrates further the importance of 
considerations other than the xl an d RMS in establishing 
the viability of a model. The xl m #4 is lower than that of 
the other POFP fits. The RMS, is as low as that of the other 
POFP cases (and lower than that of the Keplerian reproduc- 
tion, solution #2). However, when integrating over a time 
span of over 50 times that of the longer planetary period, 
the trajectory shows an ejection of one of the planets due 
to a cumulative effect of planetary interaction. This insta- 
bility is also exemplified by the negative Hill difference. The 
tradeoff between stability (in this solution, related to the 
relatively high eccentricity of the second planet) and lower 
xt requires knowledge of the former in order to evaluate the 
latter. (This aspect of the model is treated by McCarthy 
et al. (2004) when discussing the stability analysis of the 
HD160691.) 



Solution #5, with xl = 1-36 and RMS = 4.18, is an 
example of a co planar solution. Co planarity can be enforced 
during the POFP by limiting the 22 axis of the second planet 
and the tilt of its orbit, r (see appendix A). The positive Hill 
difference shows the solution to be possibly stable, which is 
in agreement with an integration over 400 times the longer 
planetary period. The agreement with the parameters for the 
first planet of solutions #1 — #2 is apparent. The period 
of the second planet in this solution is 1.56 times that of 
the second planet in solution #2. The second eccentricity in 
#5, e2 = 0.226, is about 0.4 that in #2. The combination 
of relatively low xl an d RMS, and apparent stability makes 
this a viable solution. 

Solution #6 has a xl = 1-44 and a RMS = 4.18. Unlike 
solutions #3 — #5, this one has been fitted using the entire 
14-tuple of both orbital and observational parameters (see 
§3 and appendix C). Integration of this three-dimensional 
solution over 400 times the longer planetary period shows 
an apparent stability. This is supported by the positive Hill 
difference for this case. As in the other POFP cases, the 
parameters for first planet in this solution are in agreement 
with those of McCarthy et al. (2004) found in solutions #1 — 
#2. The period of the second planet in #6 is 2.8 times that 
of the second planet in #2 while the eccentricity of this 
planet in #6 shows a nearly circular orbit. 

Solution #7, with a xl = 1-39 and RMS = 4.25, is 
not coplanar, having a relatively small angle of 5° between 
the two orbital planes. As with the previous POFP solu- 
tions, this one also agrees with the parameters published 
and reproduced, for the first planet. The second planet has 
a period 2.05 times that of the second planet in #2 and an 
eccentricity of 0.12, making for a more circular orbit. An in- 
tegration over 1000 times the longer planetary period shows 
this solution to be highly stable, a property supported by 
the positive value of the Hill difference. Fig. 1 shows the fit 
for this solution to the RV observations (where jitter has 
been included in the error bars). Fig. 2 shows a plot of the 
orbits integrated by the Solver around the CM of the sys- 
tem over a time span of 1 period of the outer planet. A plot 
of the orbit integrated over a time span of 400 times the 
period of the outer planet is shown in Fig. 3. The optimized 
look vector is shown as a red arrow stemming from the CM 
and the initial velocities as vectors stemming from the initial 
position of each orbit. 

Solutions #8 and #9, show an agreement with the pa- 
rameters for the first planet in all of the other POFP solu- 
tions. With similar periods for the second planet, they differ 
mainly in its orbital eccentricity values. The similar Hill dif- 
ference values in the two solutions attest to the comparable 
behavior in orbital stability that is shown by direct integra- 
tion of 1000 times the longer planetary period. Inspection 
of the plots of these trajectories shows the two planets in- 
teracting in a manner that is inconclusive in that it may be 
bordering instability. Using features in MATLAB that allow 
a presentation of the calculated trajectories as a function of 
time (i.e. in animation form) shows the orbits change their 
positions with respect to the star, but overall do not interact 
to the point of disruption. In that respect, it is interesting 
to note that the xl an d RMS are higher in #8, than in #9, 
possibly exemplifying higher stability which would be com- 
patible with the higher Hill difference. If stable, the nearly 
coplanar #9, with xl = 1-36 and RMS = 4.18, is a highly 
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Table 1. HD160691b (First Planet) 



Solution 


Solution 


Star Mass 


P 


c 


a 


M sin i 


MaxVr 


X. v 


RMS^ 


OV^ 


Hill # -Hill cl £ 


# 


Type 


(M ) 


(days) 




(AU) 


(Mj) 


(m/s) 




(m/s) 


(m/s) 


1 


Published^ 


1.08 


645.5 


0.200 


1.50 


1.67 


38.0 


2.56^ 
1.27/, 


4.66 






2 


Kep Rep 


1.05 


642.4 


0.200 


1.47 


1.63 


37.7 


4.20 






3 


POFP 


1.10 


642.0 


0.179 


1.53 


1.72 


36.7 


1.44~ 


4.33 


-15.4 


-0.01 


1 


POFP 


1.05 


646.7 


0.185 


1.52 


1.74 


37.8 


1.34 


4.15 


-10.2 


-0.04 


5 


POFP 


1.06 


642.6 


-0.1989 


1.49 


1.65 


37.3 


1.38 


4.18 


-22.0 


0.23 


6 


POFP 


1.08 


642.2 


0.192 _ 


1.50 


1.69 


37.7 


1.44 


4.18 


-68.9 


0.60 


7 


POFP 


1.08 


642.0 


0.183 


1.49 


1.68 


37.5 


1.39 


4.25 


-27.8 


0.43 


8 


POFP 


1.09 


643.0 


0.175 


1.50 


1.70 


37.6 


1.41 


4.27 


-21.5 


0.25 


9 


POFP 


1.08 


642.9 


0.202 


1.50 


1.70 


37.6 


1.36 


4.18 


-20.1 


0.20 



a For the generated solutions, #2 - #9, RMS = u TCS = ^/RMS| RR + cr jit 2 
" Offset Velocity (calculated as part of the xi) 

c The Hill difference which is the critical Hill stability value subtracted from the Hill stability number for 
this solution 

d McCarthy et al., 2004, ApJ 617, 575M 

e The value published is \ v = 1.6, assuming no jitter. 

/ This solution includes 3 m/s jitter. Removing the jitter results in xt = 2-7 (or Xf = 1-64). Refer to §6.1 
and §6.2. 

9 the sign of the eccentricity determines whether the apsis position (the length of which is the numerator 
in equation A3) is an apastron or periastron. 



Table 2. HD160691c (Second Planet) 



Sol' 


P 


c 


a 


Msini 


MaxVr 


h 


«2 


Comments 


# 


(days) 




(AU) 


(Mj) 


(m/s) 


(degrees) 


(degrees) 




1 


2986 


0.57 


4.17 


3.10 


51.0 






publication^ may contain an error 


2 


2977 


0.56 


4.09 


3.44 


47.7 






Double Keplerian 


3 


3616 


0.02 


4.76 


1.85 


23.0 


1.60 


52.6 


Uncertain Stability, large first-planet mass'' 


1 


5593 


0.68 


6.27 


1.75 


19.4 


1.30 


61.0 


Unstable, large first-planet mass c 


5 


4657 


0.23 


5.58 


2.73 


31.9 


12.6 


12.6 


Stable, Co-Planar 


6 


8315 


0.00 


8.26 


7.87 


74.6 


22.4 


90.0 


Stable, 3D, Optimized with 14 parameters 


7 


6131 


0.12 


6.73 


2.95 


31.0 


72.7 


76.8 


Stable, 3D (5° between orbital planes) 


8 


4991 


0.15 


5.89 


2.48 


27.7 


17.8 


69.4 


Stable (or semi-stable) 


9 


4984 


0.26 


5.88 


2.51 


28.2 


07.0 


18.4 


Stable, slightly large planetary masses 1 ^ 


















All solutions except for the publication 


















(#2 - #9) include 3 m/s jitter 



a McCarthy et al., 2004, ApJ 617, 575M 

b m pi S 62.6Mj 

c m pi 79.47Mj 

d m pi ~ 13.91 A/.,, m P2 ^ 7.95Mj 

viable solution, similar in that respect to #5 and possibly 
#7. The planetary masses in #9, are slightly larger than 
several Jupiter masses and the plane of orbit's orientation 
is nearly face on. The mass of the second planet in solution 
#8 is outside the range allowed by the KC model, as is the 
case with solution #7 (see §5.2). 



6.2 Error Analysis and Fitness-Surface 

The uncertainties for the quantities of solution #7 are shown 
in Table 6. These values were calculated by perturbing each 



of the solution's parameters and maintaining that changed 
value fixed during optimization. For each perturbation of 
a particular parameter, the search for the desired fitness 
value (x 2 not normalized) was made on the entire remain- 
ing (n — l)-tuple. Here, the size of the desired confidence 
interval between the fitness value of the original minimum 
and that of the perturbed-optimized one was taken to be 1, 
corresponding to one standard deviation in the x 2 function 
(i.e. optimizing to find Xpcrturb = Xorigmai + !)• The smallest 
perturbation for which the fitness value changes as desired 
is considered to be the uncertainty for the perturbed pa- 
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Table 3. Uncertainties for the parameters of the reproduced Kcplerian solution to 
the binary HD141929 



Model 




Parameter 


+ 


Published uncertainty 


Parameter 


(negative side 


value 


(positive side 


(for comparison) 


(and unit) 


uncertainty) 




uncertainty) 




Pi (days) 


0.02 


49.70 


0.05 


± 0.016 


Ki (ms" 1 ) 


400 


9950 


200 


± 170 


ei 


0.02 


0.393 


0.02 


± 0.008 


Lai (deg) 


5.00 


326 


2.00 


± 1.700 


Pi (days) 


0.300 


1993 


0.300 


± 0.190 



Table 4. Results for the Binary system HD141929 



Solution 


P 


To 


c 




U)2 




K 2 


RMS^ 


i 


Ki/K 2 




(days) 


(HJD-2451000) 




(deg) 


(deg) 


(m/s) 


(m/s) 


(m/s) 


(deg) 




Published^ 


49.70 


993.39 


0.393 


325.7 


145.7 


9950 


10580 


562 


~ 11 


0.94 


KC Primary^ 


49.71 


993.40 


0.403 


324.7 




9770 




518 




0.93 


KC Secondary 


49.69 


993.40 


0.412 




145.2 




10497 


365 




0.93 


POFP fit #l e 


49.69 


992.90 


0.386 


320.7 


140.7 


9508 


10119 


535 


10.74 


0.94 


POFP fit #2j_ 


49.72 


993.09 


0.406 


142.2 


322.2 


9556 


10166 


510 


10.78 


0.94 



a For all generated solutions and reproductions in this table, RMS = cr rcs = ^/RMS|, RR + Oj it 2 

^ Simultaneous fit of both stars' RV data 

c Carrier, 2002, A&A 389,475 

d individual Keplerian fit (using the KC model) 

e assuming Mi = 2.6M M 2 = 2.445M© 

f assuming Mi = 2.6M M 2 = 2.445M© 




_ 2 1 1 1 1 1 

-500 500 1000 1500 2000 2500 

time (days) 

Figure 1. The RV fit for solution #7 in Table 1. The error bars 
include a jitter of 3 m/s. 



rameter. As was the case here, the search is typically made 
using the Compass Method. Based on "The use of constant 
Chi- Square Boundaries as Confidence Limits" in Press et al. 
(1992) and on Bevington & Robinson (2003), this approach 
results in an error analysis that examines how sensitive the 
solution is to a change in each variable. It quantifies the 
shape of parameter space around a solution (Windmiller 




distance from CM (AU) 

Figure 2. The trajectories of the system in solution #7, around 
the centre of mass (the star, at the centre, is marked in red, the 
inner planet is marked in blue and the outer planet in black), 
integrated over a time span of 1 period of the outer planet. The 
red arrow at the centre is the look-vector, pointing in the direction 
of the observer's line of sight. Also shown are the initial velocities 
for each planet 

2006) and is not limited to a function resembling the one 
passing through the original data points (i.e. it reaches fam- 
ilies of the original curve). Because of that such an analysis 
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Table 5. xl intervals for different confidence levels using the F test 
with examples shown for HD160691 (u = 34) 



Confidence: 


2a 


l 

F 

la 


xl 


— > F 
la 


2a 


{ X "}min/max = Xv * 


0.493 


0.703 




1.423 


2.028 


For solution #1: 
For solution #2: 


1.310 
0.630 


1.860 
0.900 


2.650 
1.280 


3.770 
1.820 


5.370 
2.600 




Figure 3. The trajectories of the system in solution #7, into- Figure 4. The surface xl = « (ma, ma), looking down the fitness 

grated over a time span of 400 periods of the outer planet ( 703 axis > for HD160691. The fitness was calculated on a mass grid of 

years) demonstrating its stability 10 x mass values. The fitness scale is shown as a colour bar 

on the right, with a minimum value of xl min = 1-35. 



can reflect properties of the set of parameters that other- 
wise require a different way of interpretation (Examples to 
that are shown in Windmiller (2006) and Paper II). It also 
stands in contrast to methods adding or redistributing Gaus- 
sian noise, such as the Monte Carlo, frequently used in the 
literature, which, in fact, tend to underestimate the error. 
Table 4 also shows the composite errors of quantities which 
include parameters of the basic 12-tuple in their expressions. 
These were calculated in closed form using the propagation 
of errors. 

An analysis based on the translation from the published 
double Keplerian solution, into the POFP 12-tuple, was car- 
ried out by calculating the surface of fitness as a function of 
the two masses (or inclinations), for the lowest AQ in each 
grid-point (see §5.2 and Appendices C,E). Fig. 4 shows the 
surface, viewed along the xl ax i s ( m the origin's direction), 
where the different heights of the fitness are represented by 
a range of colours. Fig. 5 is a simplified contour plot of 
the main features shown in Fig 4. With the KC parame- 
ters used to generate Newtonian initial conditions for the 
Solver, the reproduced KC solution #2, with xIkc P ~ 1-27 
now has a minimum, X^Ncwt — 1-38. The xi values of the 
POFP solutions given in Table 1 are closely comparable to 
the best X^Newt values derived from the KC parameters of 
#2. Because the masses derived in a POFP solution are only 
bounded by mi,2 > 0, they are not limited to the range set 



by the KC, allowing for a solution outside the region boxed 

Student VersiorTot MATLAB 

by the fitness surface and the constant surface representing 
the single fitness value of #2. This point is discussed further 
in subsequent papers in the context of both fits for different 
stellar systems and stability analysis. 

6.3 Discussion of HD160691 Results 

Reconstructing planetary orbits by fitting RV measurements 
is an inverse problem, and as such, it is reasonable to expect 
multiple solutions that are consistent with the data. The 
variety of cases, #3 — #9, demonstrate the nature of the 
POFP operation as an attempt to solve this inverse prob- 
lem, and the types of solutions it can produce. In the inter- 
pretation presented here, the fit for every solution is judged 
by the three criteria, in order: xl> RMS, and stability. Only 
when considering all three, will the viability of a particu- 
lar solution be realized. The interchange between stability 
and lower xl i s demonstrated, for example, in the differ- 
ences between solutions #4, #5, #6, #7 increasing in stabil- 
ity with the RMS rising along with the xl (dropping again 
for #7). As discussed in §5 and §6.1, the issue of stability 
is addressed by both direct integration and the Hill sta- 
bility analysis. The Hill difference is apparent as a general 
trend but, depending on the actual value for a particular 
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Table 6. Uncertainties for the parameters of Solution #7 



Model 




Parameter 


+ 


Parameter 


(negative side 


value 


(positive side 


(and unit) 


uncertainty) 




uncertainty) 


mi (Mj) 


0.100 


1.76 


0.500 


m 2 (Mj) 


0.500 


3.03 


0.900 


m 3 (M ) 


0.080 


1.08 


0.200 


xi (AU) 


0.030 


1.22 


0.030 


x 2 (AU) 


0.300 


5.77 


0.300 


Vi (AU) 


0.300 


-1.15 


0.600 


z 2 (AU) 


0.090 


-0.48 


0.090 


6»i (rad) 


0.009 


-0.0655 


0.009 


6» 2 (rad) 


0.090 


-0.521 


0.090 


ei 


0.008 


0.183 


0.030 


e2 


0.090 


0.123 


0.090 


t (rad) 




set to t = 




Composite Error: 


mi sinii (Mj) 


0.10 


1.60 


0.50 


m2sini2 (Mj) 


0.10 


2.75 


0.50 


a x (AU) 


0.07 


1.50 


0.07 


a 2 (AU) 


0.80 


6.73 


0.80 


Pi (days) 


40 


642 


40 


P2 (days) 


500 


6130 


500 



values of their fitness function, assuming the same jitter ef- 
fect (see also Fig. 4). In any event, a full consideration of 
all three criteria seems crucial in order to assess the accept- 
ability of a model for a given solution and understand its 
implications. This full consideration is also necessary in or- 
der to compare between different solutions representing the 
same model and between different models. The observer's 
direction (expressed by the look-vector) can also be con- 
sidered. Of the six solutions suspected of being stable (or 
at least semi stable), #3, #5, #6, #7, #8, #9, solutions #5 
and #9 are nearly co-planar face on while #7 is nearly co- 
planar edge on. Solutions #3, #6, #8 are three dimensional. 
All of these solutions have eccentricity values smaller than 
0.26. The lowest fit is that of the unstable solution #4 with 
xt = 1.34. Generated for the same data set, it is interest- 
ing to note how diverse the solutions are, as the different 
observational orientations of #7 and #9 show. 

Considering the dispersion across the variety of solu- 
tions over the periods shows that the inner period is well 
constrained (and in agreement with the literature and its 
reproduction). The outer period, however, is not. Of the 
six stable-suspected solutions mentioned above, #5, #8, #9 
have similar outer periods, differing from those of #6 and 
#7 (which also differ from one another) and from the liter- 
ature and reproduction. These facts also express the inter- 
connected multitude of issues discussed in this paper, em- 
anating from the structure of POFP, especially with rela- 
tion to existing methods. Extending the time baseline of 
the observations to include sufficient data for a fit that may 
determine the period of the second planet will provide a dis- 
tinction between POFP and other models. POFP fits based 
on such data (for example the HARPS data used in Pepe 
2006) are discussed in subsequent publications. 

Examining the features of the fitness surface provides 
an insight as to the relation between the Newtonian POFP 




Figure 5. A contour representation of the main features of Fig. 
4. As shown by the height labels, the denser concentrations of 
contour lines designate lower regions of xi on the fitness surface 
(bluer regions in Fig. 4) 



solution, may not be definite. The capacity for generality 
in the geometry of the POFP solutions §"^>paxelr£ mtnis 
regard by noting that solutions #5, #7, #6 are co-planar, 
nearly co-planar and completely three dimensional respec- 
tively. Solutions #3 and #6 demonstrate the same value of 
xi, with the former being possibly unstable while the sta- 
bility of the latter being well supported. This difference in 
stability, also shown by the respective RMS values, is possi- 
bly related to the large mass of the first planet occurring in 
#3 only. Following the arguments and methods presented in 
§5, §6.1, and in this section, the disparity between a Keple- 
rian fit and a Newtonian fit is also expressed in the different 
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and the Keplerian (or KC) approximation. Of prime interest 
is the question of gravitational interaction in a physically 
realistic model. The significance of POFP in that regard, 
compared to other models, is expressed here in two ways. 
One is the variation in fitness across the (mi, 1712) domain 
of the surface, compared to the fitness value provided by the 
Keplerian solution. For example, the X^Kep value provided 
in #2 requires a specific range of masses in this domain. Al- 
ternatively, the values of m sin i provided by #2 correspond 
on the surface to higher values in X^Nowt than the x^Kep 
tabulated. The gap between the POFP and the Keplerian 
model is also exhibited with solutions existing outside the 
allowable Keplerian mass region and have a reasonable xt, 
RMS, and are stable (for example, #7 in Table 1). The high 
RCS and RMS of most of the POFP orbits (represented 
by the points of the surface) indicate that mutual planetary 
interaction plays a crucial role in considering the physics be- 
hind the models and in any comparison between them. The 
fact that viable POFP solutions, with masses not allowed 
in the Keplerian parameter space, exist, further emphasizes 
this point. 



7 FUTURE WORK 

The use of fundamental physical principles while preserving 
generality - physical, geometrical, and structural - allows the 
production of orbital solutions by POFP that are diverse, 
yet physically realistic. As demonstrated by the HDf6069f 
analysis, the use of the POFP, especially considering the va- 
riety of independent optimisation methods it is structured 
around, raises some interesting issues. Perhaps the most dis- 
tinct of those is the acceptability of the multiple (here, dou- 
ble) Keplerian model. On the one hand, comparison between 
the POFP and the Keplerian (or KC) predictions leads to 
the question of what constitutes sufficient viability of any 
model with respect to observations and also in comparison to 
another model. Considerations concerning the conditions for 
measuring a fit and the interpretation accompanying these 
conditions are highly valuable in assessing a comprehensive 
realistic description of natural phenomena. This issue has 
been expressed in this research with the requirement for 
making a wise use of the \ 2 and with arguments, such as the 
ones discussed earlier, regarding the necessity of the triple 
criteria for the acceptance of a fit (\ 2 , RMS and stability). 
These considerations demonstrate a discussion concerning 
the foundation of scientific explanations. On the other hand, 
from a physical standpoint, the analysis of POFP solutions 
or a comparison between those and their Keplerian coun- 
terparts is linked to the question of gravitational interac- 
tion between multiple orbital elements, such as planets. The 
translation of KC into Newtonian information suggests that 
the former may be incomplete, regardless of the question of 
compatibility with the observed quantities. 

In light of attempting to address those issues, at least 
two directions of utilizing the POFP will follow. One use 
is the study of other stellar systems suspected of including 
planets. An example of that is the investigation, relating to 
the existing literature, of RV data for 47 UMa and HD12661, 
which will appear in Paper II (and is also based on the exist- 
ing investigation in Windmiller (2006)). Such studies are not 
limited to Three-Body configurations as an analysis of the 



suspected multiple-planetary system p Cancri A, on Paper 
III, will show. The other avenue in which the POFP can be 
used constructively is the vast topic of stability. Assessing 
instability in each solution via a full Three-Body integration 
is computationally expensive, while short time-span integra- 
tion can only hint at such a possibility. Thus, the preliminary 
check for instability (though it may rely on founded criteria 
such as the Hill condition) like the one described for each so- 
lution here, must be expanded into a consistent methodical 
set of tools allowing for a rigorous analysis. A more thor- 
ough study of stability, based in part on the fitness surface 
described above and that may possibly point in a systematic 
direction, will be described in a future paper. 
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APPENDIX A: INITIAL POSITIONS IN THE PLANETARY MODEL 

In the orbital planes of masses mi and m.2 around mass m 3 , the initial position of each mass in orbit is: 

8ei = 1 fl ' (1 ~ ei i (cos fli, sin fli.O), s e2 = I a2(1 ~ 62 i (cos 2 , sin 02,0) (Al) 
1 + ei cos 0i 1 + e2 cos 02 

respectively. The initial velocity of each mass in orbit is: 



G{mi+mz), . . . „. G (m2 + m-i) . 

v e i = \ ^ 7- (- sm0i, ei + cos0i, 0) , v e 2 = \ j- r- (- sin 2 , e 2 + cos 2 , 0) (A2) 

V oi(l-ei) V a 2 (1 - e 2 ) 

The semi-major axes ai, 2 are given by the eccentricities ei,2 and the coordinates of the periastra, (xi) and (2:2,2/2,22) 
respectively. 



01 = -i — — , a2 = (A3) 

1 — ei 1 — e2 

The periods of orbit are given by: 



P ^ = \rT< — T (A4) 

y G (m 3 + mi, 2) 

The rotation matrices used for the transformation are: 

(cos ip — sin ip \ / cos 77 — sin ?/> \ / 1 \ 

simp cosi/> 0], R^=( 1 J , R T = ( cost -sinr J (A5) 
1 / \ sin rj cos r\ ) \ sin r cos r / 

where the angles are defined in terms of the second periastron coordinates, i.e., 
cosi/)= — X2 etc. (A6) 

The position and velocity vectors of all three bodies are then given as: 

rim = R^R^RrSei = (s e i , 0) , r 2in = R^,R„R T sf 2 , r 3in = (0,0,0) (A7) 

Wlin = R^R^RtwIi = («el, 0) , V 2 in = R^RtjRtV^, «3in = (0, 0, 0) (A8) 

The 12-tuple of orbital parameters required to generate these initial conditions in the Solver's frame is therefore C\ = 
(mi,m 2 ,m3,xi,X2,t/2,22,0i,02,ei,e 2 ,ri). 



APPENDIX B: TRANSFORMATION THROUGH THE PLANE OF THE SKY 

The Solver receives the initial conditions, Ci, specified in the previous section and produces the Three-Body trajectories 
and velocities in the form of the solution matrix xyz for time step tk 

(xi ati x~2 at 2 x~z at 3 \ 
Vl 2/1 2/2 2/2 2/3 2/3 , t k =t -* t final (Bl) 

«i ft 22 22 23 2^3 / 

(the bar represents a coordinate in the CM frame) 

One of the masses, typically m 3 , is chosen to be the star. Its velocity, v 3 = (£3,2/3,23) is extracted from the solution 
matrix. 

The observer's line of sight, or radial direction, is designated by the 2' axis which is perpendicular to the observer's plane 
of the sky. The orientation of the Solver's solution to the observer is defined by one of the orbits chosen. The angles defining 
the transformation are the orbit's longitude of periastron, w and inclination i, i.e. the 2-tuple C2 = The appropriate 

transformation matrices are the following (Hilditch 2001): 

coscj sin a; \ 

— sin uj cos lu (B2) 
1/ 



(B3) 

The velocities of equation B3 for all time steps are then compared to the observations through the fitness function of the 

x 2 . 




R* = 

and the transformation is: 
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APPENDIX C: THE LOOK VECTOR 

The solution produced by the Solver, given the orbital 12-tuple of initial conditions, Ci, is in the form of the matrix 

( (£ 3 ,y%J 3 ) \ / to \ 



V r *(t) = 



(x 3 ,y 3 ,z 3 ) ) 



corresponding to the time step vector t = 



tk 



\ (x 3 ,y 3 ,z 3 ) f J \ t f ) 

The observer-orientation from which the RV data is optimally described can be solved for. In the equation 

Vr*Li,Lj — V O 



(CI) 



/ V ol \ 

V o2 



is the vector of observed RV and Li 



cos i cos U) 
sin i sin to 
cos i 



is a 'look-vector' which designates a particular 



V V o} J 

point on a unit sphere (this point could be thought of as the direction where the line of sight leads to the observer). Because 
v r * is not a square matrix, the singular value decomposition (the MATLAB function svd) provides a convenient calculation of 
the least squares solution to this over-determined system. The svd command provides the decomposition v r „ = QSP T . With 
n being the number of rows in v r * , Q is an n x n orthogonal matrix, P is a 3 x 3 orthogonal matrix and S is a diagonal n x 3 
matrix. The matrix S _1 is constructed as a matrix of zeros, with three diagonal elements which are the reciprocals of the 
elements in S. The look- vector is solved for as 

L ilU = PS 1 Q T Vo 



From the look-vector, i and u> (i.e. the components of the directional initial conditions 2-tuple, the vector C2) are obtained. 



APPENDIX D: X AS A FITNESS FUNCTION 



Given the observed velocities (Vo)k, their errors &k, the computed (predicted) velocities {Vc)k and the offset velocities (Vs)k 
along the indices k = ki, ki ■ ■ ■ that divide the data set, the form used for the fitness function (x 2 ) is: 



ffit (Vo k ,(Tk,Vc k ,Vs k ) 



E 

ki L 



Vo kl - (Vc kl +Vs kl ] 



E 



Vo k2 - (V Ck2 + V Sk2 



+ 



(Dl) 



The indices are determined by the number of different instruments used to obtain the observations. The total number of 
entries (k, or the number of indexed points fci + k-2 + . . .), on which the summation in Equation Dl is made, is the number 
of data points. To optimize for the offset velocity, the partial derivatives of the fitness function are taken and set to zero: 



k 



9V Sk 



Vo k - (V Ck + V Sk 



Ok 







E 



Vo k - V Ck 



arriving at the expression for the optimal offset correction 

~v 0k -v c - 



Vs k = 



(D2) 



APPENDIX E: TRANSLATING A DOUBLE KEPLERIAN INTO POFP 

In the following process, the standard notation for the orbital properties of binaries, expressed in Hilditch (2001), is used. 
The orbital-observational parameters characterizing an orbit in the Keplerian model are P, K,T,u>,e. 

Given the times of first observation, to 1 , to 2 , the times of periastron passage for each planet, T\,T^ can be translated 
into the delay angles Q\ , 62 respectively. This can be done by using the Newton- Raphson method to iterate for the eccentric 
anomaly, E, in Kepler's Equation (3): 

2 ^ toi, 2 1£ _ E ^ ^ - ei 2 s i n E 0l 2 (El) 

The eccentric anomaly can then be used to solve for the true anomaly at time to (4): 
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di,2 = 0(ioi, 2 ) = 2arctan 



1 + ei, 
1-ei, 



■ tan 



E 0l<2 



(E2) 



The periastron position in the first-planet coordinate system is set on the positive x axis and therefore can be calculated from 
the period (A4): 



ai,2 



PifiG (m.3 + rrn, 2) 



4tt 2 



01 (l - el) 

xi = 01 (1 - ei) = — H-77^ cos6, i(0) 

1 + eicost/i(0) 



(E3) 



where the delay angle is 6 l i(0) = 0. The two masses are either given directly or can be calculated from the inclinations of the 
two orbits by use of the mass function (which is given and by the fitted half-amplitude velocity K, period P and eccentricity 
e): 



f{m)i,2 = /(mi, 2 ) = \A ~ e i,2 > /( m )i,2 (mi, 2 + m a ) 2 = m?, 2 sin 3 ii, 2 

where mi, 2 can be extracted from the expression on the right using a Newton-Raphson iteration. 



(E4) 



The coordinates of the second periastron can be calculated using the relative positions of the two orbits' periastra. The 
rotation matrices used to transform the coordinate system of the first periastron to that of the second are: 



/ 1 \ 
Rj = I cos i — sin i , Rq - 
\ sini cosi / 

with f2i,2 being the position angle of the nodal point in each orbit. 
The transformation is then: 




cos O — sin 57 
sin Q cos 57 
1 



(E5) 



X2 
V2 

Z2 



— R— R— ii R— sii Rn 2 "»2 



02(1 - e 2 ) 





— R_ W1 R_ij R(n 2 _n 1 )Ri 2 R W2 



02(1 - e 2 ) 





(E6) 



Finally, the tilt r, of the second orbit, relative to the first, is calculated via the same transformation, using the third 
component of the vector: 



r — R- L j 1 R-tj R(n 2 _n 1 ) Ri 2 Ru> 2 




arctan (r z ) 



The look-vector can also be calculated with a transformation using the same rotation matrices: 



(E7) 



(E8) 



This paper has been typeset from a Tj^X/ ETp^X file prepared by the author. 



